### Load required library ###
library(lme4)
library(texreg)
library(effects)
library(openintro)

### Load data ###
# Download data from Harvard dataverse
# Set working directory, change line below as as appropriate
setwd("~/Library/CloudStorage/Dropbox/Research/WTO paper/Data")

# Load data
load("data_what_drives_international_cooperation.Rdata")


### Codebook ###

# Two datasets will appear in the global environment:

# wto: dyadic and annual data, used for all models in Table 1 and 2 with the exception of the last model of Table 1
# wto2: dyadic data with merged years, used for the last model of Table 1

# The followng variables are found in both datasets, in the wto dataset they are annual, in the wto2 dataset combined over the years. For a more detailed description, see methods section of the paper:

# iso3c1: ISO3 code of first country in dyad
# iso3c2: ISO3 code of second country in dyad
# CoopTotalCorr: cooperative events between dyads in question (by year for wto dataset, over all years for wto2 dataset)
# CoopAG: cooperative events for agricultural issue only (only in wto dataset)
# CoopNAMA: cooperative events for NAMA issue only (only in wto dataset)
# gwf: dyadic regime combination (10 categories, see paper)
# gwf2: dyadic regime combinations (only 3 categories: two democracies, mixed dyads, two autocracies)
# logtrade: bilateral trade flows between dyad in question, logged and lagged
# logdist: geographic proximity between dyad in question, logged
# diffGDP: difference between the GDP of the two countries in question, logged and lagged
# totalGDP: combined GDP of the two countries in question, logged and lagged
# comlang_off: cultural proximity measued by common official language
# groups: is dyad in question in one (or more) joind negotiation groups
# Year: year of observation


### Notes for models:
# Selection stage models are slow and will take a while to run
# Year dummies in models are omitted from the paper

### Models Table 1 ###
# Note: the year effects are excluded from the tables in the paper

# Selection stage model
mod1_sel <- glmer(CoopTotalCorr>0 ~ gwf2 * logtrade + logdist + diffGDP + totalGDP + comlang_off + colony + groups + Year + (1 | iso3c2) + (1 | iso3c1), data=wto, family=binomial(link = "logit"))
summary(mod1_sel)

# Allocation stage model 
mod1_alloc <- lmer(CoopTotalCorr ~ gwf2 * logtrade + logdist + diffGDP + totalGDP + comlang_off + colony +  groups + Year + (1 | iso3c2) + (1 | iso3c1), data=wto, subset =  CoopTotalCorr>0)
summary(mod1_alloc)

# Full model 
mod1_full <- lmer(CoopTotalCorr ~ gwf2 * logtrade + logdist + diffGDP + totalGDP + comlang_off + colony + groups + Year + (1 | iso3c2) + (1 | iso3c1), data=wto)
summary(mod1_full)

# Full model (years combined)
mod1_merge <- lmer(CoopTotalCorr ~ gwf2 * logtrade + logdist  + diffGDP + totalGDP + comlang_off + colony +  groups + (1 | iso3c2) + (1 | iso3c1), data=wto2)
summary(mod1_merge)


# Make Table 2
screenreg(list(mod1_sel, mod1_alloc, mod1_full, mod1_merge), digits = 3)


### Models Table 2 ###
# Note: the year effects are excluded from the tables in the paper

# Selection stage model 
mod2_sel <- glmer(CoopTotalCorr>0 ~ gwf * logtrade + logdist + diffGDP + totalGDP + comlang_off + colony + groups + Year + (1 | iso3c2) + (1 | iso3c1), data=wto, family=binomial(link = "logit"))
summary(mod2_sel)

# Allocation stage 
mod2_alloc <- lmer(CoopTotalCorr ~ gwf * logtrade + logdist + diffGDP + totalGDP + comlang_off + colony +  groups + Year + (1 | iso3c2) + (1 | iso3c1), data=wto, subset =  CoopTotalCorr>0)
summary(mod2_alloc)

# Full model 
mod2_full <- lmer(CoopTotalCorr ~ gwf * logtrade + logdist + diffGDP + totalGDP + comlang_off + colony + groups + Year + (1 | iso3c2) + (1 | iso3c1), data=wto)
summary(mod2_full)

# Full model (years combined)
mod2_merge <- lmer(CoopTotalCorr ~ gwf * logtrade + logdist  + diffGDP + totalGDP + comlang_off + colony +  groups + (1 | iso3c2) + (1 | iso3c1), data=wto2)
summary(mod2_merge)

# Agriculture selection model
mod_ag_sel <- glmer(CoopAg>0 ~ gwf * logtrade + logdist + diffGDP + totalGDP + comlang_off + colony + groups + Year  + (1 | iso3c2) + (1 | iso3c1), data=wto, family=binomial(link = "logit"))
summary(mod_ag_sel)

# Agriculture allocation model
mod_ag_alloc <- lmer(CoopAg ~ gwf * logtrade + logdist + diffGDP + totalGDP + comlang_off + colony + groups + Year  + (1 | iso3c2) + (1 | iso3c1), data=wto,  subset = CoopAg>0)
summary(mod_ag_alloc)

# NAMA selection model
mod_nama_sel <- glmer(CoopNAMA>0 ~ gwf * logtrade + logdist + diffGDP + totalGDP + comlang_off + colony + groups + Year  + (1 | iso3c2) + (1 | iso3c1), data=wto, family=binomial(link = "logit"))
summary(mod_nama_sel)

# NAMA allocation model
mod_nama_alloc <- lmer(CoopNAMA ~ gwf * logtrade + logdist + diffGDP + totalGDP + comlang_off + colony + groups + Year  + (1 | iso3c2) + (1 | iso3c1), data=wto,  subset = CoopNAMA>0)
summary(mod_nama_alloc)

# Make Table 2
screenreg(list(mod2_sel, mod2_alloc, mod2_full, mod2_merge, mod_ag_sel, mod_ag_alloc, mod_nama_sel, mod_nama_alloc), digits = 3)



### Figure 1 ###
# Code will produce plot in pdf and write it into working directory

myPDF("Figure1.pdf",10,5.5)
par(mfrow=c(1,2), mar=c(3,3,2,1),mgp=c(1.3,.5,0),cex=1.1, las=0)

# Selection stage figure
x <- data.frame(effect(c("gwf2", "logtrade"), mod1_sel,confidence.level=.95, xlevels=list(logtrade=seq(0,13.5,.25))))


y <- x[x$gwf2=="2dems",]
plot(y$fit ~ y$logtrade, type="l", ylim=c(0,1),
     ylab="Predicted prob. of cooperation", xlab="Total trade (in billion US$)", xaxt='n')
axis(1, at=log(c(0,1,10,100,1000,10000,100000)+1), labels = c(0,1,10,100,1000,10000,"100000"))
axis(1, at=log(c(1,100000)+1), labels = c(1,"100000"))
abline(h=seq(0,1,.2), lty=2, lwd=1, col="lightgray")
abline(v=log(c(0,1,10,100,1000,10000,100000,500000)+1), lty=2, lwd=1, col="lightgray")

polygon(c(y$logtrade, rev(y$logtrade)), c(y$lower, rev(y$upper)), col="#808080D6", border = NA)

y <- x[x$gwf2=="mixed" & x$logtrade<13.45,]
polygon(c(y$logtrade, rev(y$logtrade)), c(y$lower, rev(y$upper)), col="#C3C3C399", border = NA)

y <- x[x$gwf=="2aut" & x$logtrade<11.5,]
polygon(c(y$logtrade, rev(y$logtrade)), c(y$lower, rev(y$upper)), col="#A9A9A999", border = NA)

y <- x[x$gwf=="2dems",]
lines(y$fit ~ y$logtrade, type="l", lwd=2)

y <- x[x$gwf=="mixed" & x$logtrade<13.5,]
lines(y$fit ~ y$logtrade, type="l", lwd=2, lty=2)

y <- x[x$gwf=="2aut" & x$logtrade<11.5,]
lines(y$fit ~ y$logtrade, type="l", lwd=2, lty=3)

legend("topleft", c("2 Democracies", "Mixed dyads", "2 Autocracies"), lty=c(1,2,3), bty="n", cex = .8, lwd = 2)
title("a) Selection Stage", adj=0)

# Allocation stage figure
x <- data.frame(effect(c("gwf2", "logtrade"), mod1_alloc,confidence.level=.9, xlevels=list(logtrade=seq(0,13.5,.25))))

y <- x[x$gwf2=="2dems",]
plot(y$fit ~ y$logtrade, type="l", ylim=c(1,5),
     ylab="Predicted prob. of cooperation", xlab="Total trade (in billion US$)", xaxt='n')
axis(1, at=log(c(0,1,10,100,1000,10000,100000)+1), labels = c(0,1,10,100,1000,10000,"100000"))
axis(1, at=log(c(1,100000)+1), labels = c(1,"100000"))
abline(h=seq(0,5,1), lty=2, lwd=1, col="lightgray")
abline(v=log(c(0,1,10,100,1000,10000,100000,500000)+1), lty=2, lwd=1, col="lightgray")

polygon(c(y$logtrade, rev(y$logtrade)), c(y$lower, rev(y$upper)), col="#808080D6", border = NA)

y <- x[x$gwf2=="mixed" & x$logtrade<13.45,]
polygon(c(y$logtrade, rev(y$logtrade)), c(y$lower, rev(y$upper)), col="#C3C3C399", border = NA)

y <- x[x$gwf=="2aut" & x$logtrade<11.5,]
polygon(c(y$logtrade, rev(y$logtrade)), c(y$lower, rev(y$upper)), col="#A9A9A999", border = NA)

y <- x[x$gwf=="2dems",]
lines(y$fit ~ y$logtrade, type="l", lwd=2)

y <- x[x$gwf=="mixed" & x$logtrade<13.5,]
lines(y$fit ~ y$logtrade, type="l", lwd=2, lty=2)

y <- x[x$gwf=="2aut" & x$logtrade<11.5,]
lines(y$fit ~ y$logtrade, type="l", lwd=2, lty=3)

legend("topleft", c("2 Democracies", "Mixed dyads", "2 Autocracies"), lty=c(1,2,3), bty="n", cex = .8, lwd = 2)
title("a) Allocation Stage", adj=0)
dev.off()




### Figure 2 ###
# Code will produce plot in pdf and write it into working directory

myPDF("Figure2.pdf",15,5.5)
par(mfrow=c(1,3), mar=c(3,3,2,1),mgp=c(1.3,.5,0),cex=1.1, las=0)
x <- data.frame(effect("gwf * logtrade",mod2_sel,confidence.level=.9, xlevels=list(logtrade=seq(0,13.5,.25))))

# 4 equal dyads

y <- x[x$gwf=="2dems",]
plot(y$fit ~ y$logtrade, type="l", ylim=c(0,1),
     ylab="Predicted number of cooperative ties", xlab="Total trade (in billion US$)", xaxt='n')
axis(1, at=log(c(0,1,10,100,1000,10000,100000)+1), labels = c(0,1,10,100,1000,10000,"100000"))
axis(1, at=log(c(1,100000)+1), labels = c(1,"100000"))
abline(h=seq(0,1,.2), lty=2, lwd=1, col="lightgray")
abline(v=log(c(0,1,10,100,1000,10000,100000,500000)+1), lty=2, lwd=1, col="lightgray")

polygon(c(y$logtrade, rev(y$logtrade)), c(y$lower, rev(y$upper)), col="#808080D6", border = NA)

y <- x[x$gwf=="2per" & x$logtrade<8.5,]
polygon(c(y$logtrade, rev(y$logtrade)), c(y$lower, rev(y$upper)), col="#C3C3C399", border = NA)

y <- x[x$gwf=="2par" & x$logtrade<11.5,]
polygon(c(y$logtrade, rev(y$logtrade)), c(y$lower, rev(y$upper)), col="#A9A9A999", border = NA)

y <- x[x$gwf=="2oth" & x$logtrade<9.5,]
polygon(c(y$logtrade, rev(y$logtrade)), c(y$lower, rev(y$upper)), col="#E7E7E799", border = NA)

y <- x[x$gwf=="2dems",]
lines(y$fit ~ y$logtrade, type="l", lwd=2)

y <- x[x$gwf=="2per" & x$logtrade<8.5,]
lines(y$fit ~ y$logtrade, type="l", lwd=2, lty=2)

y <- x[x$gwf=="2par" & x$logtrade<11.5,]
lines(y$fit ~ y$logtrade, type="l", lwd=2, lty=3)

y <- x[x$gwf=="2oth" & x$logtrade<9.5,]
lines(y$fit ~ y$logtrade, type="l", lwd=2, lty=4)

legend("topleft", c("2 Democracies", "2 Personalist Regimes", "2 Party Regimes" ,"2 Military Regime"), lty=c(1,2,3,4), bty="n", cex = .8, lwd = 2)
title("a) Same regime dyads", adj=0)


# 3 dyads with dem

y <- x[x$gwf=="dem-per"& x$logtrade<11.1,]
plot(y$fit ~ y$logtrade, type="l", ylim=c(0,1), xlim=c(0,13.5),
     ylab="Predicted number of cooperative ties", xlab="Total trade (in billion US$)", xaxt='n')
axis(1, at=log(c(0,1,10,100,1000,10000,100000)+1), labels = c(0,1,10,100,1000,10000,"100000"))
axis(1, at=log(c(1,100000)+1), labels = c(1,"100000"))
abline(h=seq(0,1,.2), lty=2, lwd=1, col="lightgray")
abline(v=log(c(0,1,10,100,1000,10000,100000,500000)+1), lty=2, lwd=1, col="lightgray")

polygon(c(y$logtrade, rev(y$logtrade)), c(y$lower, rev(y$upper)), col="#808080D6", border = NA)

y <- x[x$gwf=="dem-par",]
polygon(c(y$logtrade, rev(y$logtrade)), c(y$lower, rev(y$upper)), col="#C3C3C399", border = NA)

y <- x[x$gwf=="dem-oth" & x$logtrade<11.3,]
polygon(c(y$logtrade, rev(y$logtrade)), c(y$lower, rev(y$upper)), col="#E7E7E799", border = NA)

y <- x[x$gwf=="dem-per" & x$logtrade<11.1,]
lines(y$fit ~ y$logtrade, type="l", lwd=2)

y <- x[x$gwf=="dem-par",]
lines(y$fit ~ y$logtrade, type="l", lwd=2, lty=2)

y <- x[x$gwf=="dem-oth" & x$logtrade<11.3,]
lines(y$fit ~ y$logtrade, type="l", lwd=2, lty=3)

legend("topleft", c("Democracy - Personalist Regimes", "Democracy - Party Regimes", "Democracy - Military Regimes"), lty=c(1,2,3), bty="n", cex = .8, lwd = 2)
title("b) Mixed dyads (one democracy)", adj=0)

# 3 dyads without dems

y <- x[x$gwf=="par-per"& x$logtrade<9.6,]
plot(y$fit ~ y$logtrade, type="l", ylim=c(0,1), xlim=c(0,13.5),
     ylab="Predicted number of cooperative ties", xlab="Total trade (in billion US$)", xaxt='n')
axis(1, at=log(c(0,1,10,100,1000,10000,100000)+1), labels = c(0,1,10,100,1000,10000,"100000"))
axis(1, at=log(c(1,100000)+1), labels = c(1,"100000"))
abline(h=seq(0,1,.2), lty=2, lwd=1, col="lightgray")
abline(v=log(c(0,1,10,100,1000,10000,100000,500000)+1), lty=2, lwd=1, col="lightgray")
title("c) Mixed dyads (no democracy)", adj=0)

polygon(c(y$logtrade, rev(y$logtrade)), c(y$lower, rev(y$upper)), col="#808080D6", border = NA)

y <- x[x$gwf=="par-oth" & x$logtrade<11.1,]
polygon(c(y$logtrade, rev(y$logtrade)), c(y$lower, rev(y$upper)), col="#C3C3C399", border = NA)

y <- x[x$gwf=="per-oth" & x$logtrade<6.7,]
polygon(c(y$logtrade, rev(y$logtrade)), c(y$lower, rev(y$upper)), col="#E7E7E799", border = NA)

y <- x[x$gwf=="par-per" & x$logtrade<9.6, ]
lines(y$fit ~ y$logtrade, type="l", lwd=2)

y <- x[x$gwf=="par-oth" & x$logtrade<11.1, ]
lines(y$fit ~ y$logtrade, type="l", lwd=2, lty=2)

y <- x[x$gwf=="per-oth" & x$logtrade<6.7,]
lines(y$fit ~ y$logtrade, type="l", lwd=2, lty=3)

#axis(1, at=log(c(0,1,10,100,1000,10000,100000)), labels = c(-1,0,10,100,1000,10000,"100000"))
legend("topleft", c("Party - Personalist Regimes", "Party - Military Regimes", "Personalist - Military Regimes"), lty=c(1,2,3), bty="n", cex = .8, lwd = 2)
dev.off()




### Figure 3 ###
# Code will produce plot in pdf and write it into working directory

myPDF("Figure3.pdf",15,5.5)
par(mfrow=c(1,3), mar=c(3,3,2,1),mgp=c(1.3,.5,0),cex=1.1, las=0)
x <- data.frame(effect("gwf * logtrade",mod2_alloc ,confidence.level=.9, xlevels=list(logtrade=seq(0,13.5,.25))))

# 4 equal dyads

y <- x[x$gwf=="2dems",]
plot(y$fit ~ y$logtrade, type="l", ylim=c(0,5),
     ylab="Predicted number of cooperative ties", xlab="Total trade (in billion US$)", xaxt='n')
axis(1, at=log(c(0,1,10,100,1000,10000,100000)+1), labels = c(0,1,10,100,1000,10000,"100000"))
axis(1, at=log(c(1,100000)+1), labels = c(1,"100000"))
abline(h=seq(0,5,1), lty=2, lwd=1, col="lightgray")
abline(v=log(c(0,1,10,100,1000,10000,100000,500000)+1), lty=2, lwd=1, col="lightgray")

polygon(c(y$logtrade, rev(y$logtrade)), c(y$lower, rev(y$upper)), col="#808080D6", border = NA)

y <- x[x$gwf=="2per" & x$logtrade<8.5,]
polygon(c(y$logtrade, rev(y$logtrade)), c(y$lower, rev(y$upper)), col="#C3C3C399", border = NA)

y <- x[x$gwf=="2par" & x$logtrade<11.5,]
polygon(c(y$logtrade, rev(y$logtrade)), c(y$lower, rev(y$upper)), col="#A9A9A999", border = NA)

y <- x[x$gwf=="2oth" & x$logtrade<9.5,]
polygon(c(y$logtrade, rev(y$logtrade)), c(y$lower, rev(y$upper)), col="#E7E7E799", border = NA)

y <- x[x$gwf=="2dems",]
lines(y$fit ~ y$logtrade, type="l", lwd=2)

y <- x[x$gwf=="2per" & x$logtrade<8.5,]
lines(y$fit ~ y$logtrade, type="l", lwd=2, lty=2)

y <- x[x$gwf=="2par" & x$logtrade<11.5,]
lines(y$fit ~ y$logtrade, type="l", lwd=2, lty=3)

y <- x[x$gwf=="2oth" & x$logtrade<9.5,]
lines(y$fit ~ y$logtrade, type="l", lwd=2, lty=4)

legend("topleft", c("2 Democracies", "2 Personalist Regimes", "2 Party Regimes" ,"2 Military Regime"), lty=c(1,2,3,4), bty="n", cex = .8, lwd = 2)
title("a) Same regime dyads", adj=0)


# 3 dyads with dem

y <- x[x$gwf=="dem-per"& x$logtrade<11.1,]
plot(y$fit ~ y$logtrade, type="l", ylim=c(0,5), xlim=c(0,13.5),
     ylab="Predicted number of cooperative ties", xlab="Total trade (in billion US$)", xaxt='n')
axis(1, at=log(c(0,1,10,100,1000,10000,100000)+1), labels = c(0,1,10,100,1000,10000,"100000"))
axis(1, at=log(c(1,100000)+1), labels = c(1,"100000"))
abline(h=seq(0,5,1), lty=2, lwd=1, col="lightgray")
abline(v=log(c(0,1,10,100,1000,10000,100000,500000)+1), lty=2, lwd=1, col="lightgray")

polygon(c(y$logtrade, rev(y$logtrade)), c(y$lower, rev(y$upper)), col="#808080D6", border = NA)

y <- x[x$gwf=="dem-par",]
polygon(c(y$logtrade, rev(y$logtrade)), c(y$lower, rev(y$upper)), col="#C3C3C399", border = NA)

y <- x[x$gwf=="dem-oth" & x$logtrade<11.3,]
polygon(c(y$logtrade, rev(y$logtrade)), c(y$lower, rev(y$upper)), col="#E7E7E799", border = NA)

y <- x[x$gwf=="dem-per" & x$logtrade<11.1,]
lines(y$fit ~ y$logtrade, type="l", lwd=2)

y <- x[x$gwf=="dem-par",]
lines(y$fit ~ y$logtrade, type="l", lwd=2, lty=2)

y <- x[x$gwf=="dem-oth" & x$logtrade<11.3,]
lines(y$fit ~ y$logtrade, type="l", lwd=2, lty=3)

legend("topleft", c("Democracy - Personalist Regimes", "Democracy - Party Regimes", "Democracy - Military Regimes"), lty=c(1,2,3), bty="n", cex = .8, lwd = 2)
title("b) Mixed dyads (one democracy)", adj=0)

# 3 dyads without dems

y <- x[x$gwf=="par-per"& x$logtrade<9.6,]
plot(y$fit ~ y$logtrade, type="l", ylim=c(0,5), xlim=c(0,13.5),
     ylab="Predicted number of cooperative ties", xlab="Total trade (in billion US$)", xaxt='n')
axis(1, at=log(c(0,1,10,100,1000,10000,100000)+1), labels = c(0,1,10,100,1000,10000,"100000"))
axis(1, at=log(c(1,100000)+1), labels = c(1,"100000"))
abline(h=seq(0,5,1), lty=2, lwd=1, col="lightgray")
abline(v=log(c(0,1,10,100,1000,10000,100000,500000)+1), lty=2, lwd=1, col="lightgray")
title("c) Mixed dyads (no democracy)", adj=0)

polygon(c(y$logtrade, rev(y$logtrade)), c(y$lower, rev(y$upper)), col="#808080D6", border = NA)

y <- x[x$gwf=="par-oth" & x$logtrade<11.1,]
polygon(c(y$logtrade, rev(y$logtrade)), c(y$lower, rev(y$upper)), col="#C3C3C399", border = NA)

y <- x[x$gwf=="per-oth" & x$logtrade<6.7,]
polygon(c(y$logtrade, rev(y$logtrade)), c(y$lower, rev(y$upper)), col="#E7E7E799", border = NA)

y <- x[x$gwf=="par-per" & x$logtrade<9.6, ]
lines(y$fit ~ y$logtrade, type="l", lwd=2)

y <- x[x$gwf=="par-oth" & x$logtrade<11.1, ]
lines(y$fit ~ y$logtrade, type="l", lwd=2, lty=2)

y <- x[x$gwf=="per-oth" & x$logtrade<6.7,]
lines(y$fit ~ y$logtrade, type="l", lwd=2, lty=3)

#axis(1, at=log(c(0,1,10,100,1000,10000,100000)), labels = c(-1,0,10,100,1000,10000,"100000"))
legend("topleft", c("Party - Personalist Regimes", "Party - Military Regimes", "Personalist - Military Regimes"), lty=c(1,2,3), bty="n", cex = .8, lwd = 2)
dev.off()

